function [results] = XPTmeter_MASTER(data, TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,garnet_coefficients, spinel_coefficients, plagioclase_coefficients)
%fits,dataX,garnet_modelX,spinel_modelX,plag_modelX
% Elements,Weights,VariablesT,VariablesX,VariablesP,

%The process is to:
%(1) first calculate the pressure
%(2) then using that P calculate T and X
%(3) calculate the X difference to see if the liquid is is EQ with a
%lherzolite

%%%%%%%%%%%%%

% [A,TemperatureIndicies] = ismember(VariablesT, Elements);
% [A,CompIndicies] = ismember(VariablesX, Elements);
% [A,PressureIndicies] = ismember(VariablesP, Elements);
% 
% 
 FittedTargetVariables = {'Qtz' 'Plag' 'Oliv' 'Cpx'};
% [A,FittedVariablesIndicies] = ismember(FittedTargetVariables, Elements);
% 
% 
 ModelVariables = {'T' 'Oliv' 'Cpx' 'Plag' 'Qtz'};
% [A,ModelVariablesIndicies] = ismember(FittedTargetVariables, ModelVariables);
% 


data(isnan(data))=0;
onesmatrix = ones(size(data,1),1);

if nargin > 1
    
    %rows are:
    % T
    % oliv
    % cpx
    % plag
    % qtz
    % P
    
    
    %Step 1: calculate P
    variableValues = [onesmatrix data(:,PressureIndicies)]; %add in a ones column account for the intercept
    %variableValuesGarnet = [onesmatrix data(:,PressureIndiciesGarnet)]; %add in a ones column account for the intercept
    
    garnetPressures = variableValues*garnet_coefficients(end,:)'; %in kbar
    spinelPressures = variableValues*spinel_coefficients(end,:)'; %in kbar
    plagioclasePressures = variableValues*plagioclase_coefficients(end,:)';%in kbar
    
    
    %Step 2 : calculate T and X
    
    
    variableValues = [onesmatrix (garnetPressures) data(:,TemperatureIndicies)];
    T = [variableValues*garnet_coefficients(1:end-1,:)'];
    variableValues = [onesmatrix (garnetPressures) data(:,CompIndicies)];
    X = [variableValues*garnet_coefficients(1:end-1,:)'];
    garnetTX=[T(:,1), X(:,2:end)];
    
    
    variableValues = [onesmatrix spinelPressures data(:,TemperatureIndicies)];
    T = [variableValues*spinel_coefficients(1:end-1,:)'];
    variableValues = [onesmatrix (spinelPressures) data(:,CompIndicies)];
    X = [variableValues*spinel_coefficients(1:end-1,:)'];
    spinelTX=[T(:,1), X(:,2:end)];
    
    
    variableValues =  [onesmatrix plagioclasePressures data(:,TemperatureIndicies)]; %[onesmatrix (plagioclasePressures)' repmat(data(i,TemperatureIndicies),size(onesmatrix,1),1)];
    T = [variableValues*plagioclase_coefficients(1:end-1,:)'];
    variableValues = [onesmatrix (plagioclasePressures) data(:,CompIndicies)];
    X = [variableValues*plagioclase_coefficients(1:end-1,:)'];
    plagioclaseTX=[T(:,1), X(:,2:end)];
    
    
    %%
    %Step 3 : compare the data to the differet model melt components
    
    %NOTE: tried using using the AIT distance, but doesn't work becuase Tormey components can be negative
    
    dataX = data(:,FittedVariablesIndicies);
    garnet_modelX = garnetTX(:,ModelVariablesIndicies);
    spinel_modelX = spinelTX(:,ModelVariablesIndicies);
    plag_modelX = plagioclaseTX(:,ModelVariablesIndicies);

   
    dataX_norm=bsxfun(@rdivide,dataX,nansum(dataX,2));
    garnet_modelX_norm = bsxfun(@rdivide,garnet_modelX,nansum(garnet_modelX,2));
    spinel_modelX_norm = bsxfun(@rdivide,spinel_modelX,nansum(spinel_modelX,2));
    plag_modelX_norm = bsxfun(@rdivide,plag_modelX,nansum(plag_modelX,2));
    
    
    %the same old same old way of fitting
    garnetRMSD=       sqrt(sum((garnet_modelX_norm-dataX_norm).^2, 2)./size(FittedTargetVariables,2));
    
 
    spinelRMSD=       sqrt(sum((spinel_modelX_norm-dataX_norm).^2, 2)./size(FittedTargetVariables,2));
    plagioclaseRMSD=  sqrt(sum((plag_modelX_norm-dataX_norm).^2, 2)./size(FittedTargetVariables,2));
    
    
    fits=[garnetRMSD spinelRMSD plagioclaseRMSD];
    

    %%

    results.dataX = dataX;   %the melt components of the data
    results.garnet_modelX = garnet_modelX;  %the melt components of the garnet model
    results.spinel_modelX = spinel_modelX; %the melt components of the garnet model
    results.plag_modelX = plag_modelX; %the melt components of the garnet model
    
    
    %Now calculate fits for all combinations of 3:
%     n=1:4;
%     c=nchoosek(n,3);
%     for d=1:size(c,1)
%         
%         componentsHere = c(d,:);
%         FittedTargetVariables_sub = FittedTargetVariables(componentsHere);
%         
%         [A,FittedVariablesIndicies] = ismember(FittedTargetVariables_sub, Elements);
%         [A,ModelVariablesIndicies] = ismember(FittedTargetVariables_sub, ModelVariables);
%         
%         dataX = data(:,FittedVariablesIndicies);
%         garnet_modelX = garnetTX(:,ModelVariablesIndicies);
%         spinel_modelX = spinelTX(:,ModelVariablesIndicies);
%         plag_modelX = plagioclaseTX(:,ModelVariablesIndicies);
%         
%         dataX_norm=bsxfun(@rdivide,dataX,nansum(dataX,2));
%         garnet_modelX_norm = bsxfun(@rdivide,garnet_modelX,nansum(garnet_modelX,2));
%         spinel_modelX_norm = bsxfun(@rdivide,spinel_modelX,nansum(spinel_modelX,2));
%         plag_modelX_norm = bsxfun(@rdivide,plag_modelX,nansum(plag_modelX,2));
%           
%         %the same old same old way of fitting
%         garnetRMSD=       sqrt(sum((garnet_modelX_norm-dataX_norm).^2, 2)./size(FittedTargetVariables,2));
%         spinelRMSD=       sqrt(sum((spinel_modelX_norm-dataX_norm).^2, 2)./size(FittedTargetVariables,2));
%         plagioclaseRMSD=  sqrt(sum((plag_modelX_norm-dataX_norm).^2, 2)./size(FittedTargetVariables,2));
%         
%         fits_5(:,d,:) = [garnetRMSD spinelRMSD  plagioclaseRMSD];     
%     end
%      
%     
%     Weights=Weights./sum(Weights);
%     
%     %     {'Qtz'}    {'Plag'}    {'Oliv'}
%     %     {'Qtz'}    {'Plag'}    {'Cpx'} dont really use
%     %     {'Qtz'}    {'Oliv'}    {'Cpx'}
%     %     {'Plag'}    {'Oliv'}    {'Cpx'}
%     
%     for hmn = 1:size(fits_5,1)
%         AllFits = [fits(hmn,:); squeeze(fits_5(hmn,:,:))];
%         Fits_5_combined(hmn,:) = (AllFits'*Weights')';
%     end

    %%
    % modelX_normalized = bsxfun(@rdivide,modelX,nansum(modelX,2));
    % dataX_normalized = bsxfun(@rdivide,dataX,nansum(dataX,2));
    % NRMSD= sqrt(sum((modelX_normalized-dataX_normalized).^2, 2));
    %  fits = [RMSD,NRMSD] ;
    % figure(); hold on; plot(NRMSD,RMSD,'*k')
    % xlabel('NRSMD')
    % ylabel('RMSD')
    
    %%

    %results.fits = Fits_5_combined;  %contains the fit values in order garnet, spinel, plag
    results.fits = fits;  %contains the fit values in order garnet, spinel, plag
    
    adiabat = 1.5;
    
%     results.PTTmp = [...
%         garnetPressures garnetTX(:,1) garnetTX(:,1) - adiabat.*garnetPressures...
%         spinelPressures_kbars spinelTX(:,1) spinelTX(:,1) - adiabat.*spinelPressures_kbars...
%         plagioclasePressures_kbars plagioclaseTX(:,1) plagioclaseTX(:,1) - adiabat.*plagioclasePressures_kbars];
%     
     results.PTTmp = [...
        garnetPressures garnetTX(:,1) garnetTX(:,1) - adiabat.*garnetPressures...
        spinelPressures spinelTX(:,1) spinelTX(:,1) - adiabat.*spinelPressures...
        plagioclasePressures plagioclaseTX(:,1) plagioclaseTX(:,1) - adiabat.*plagioclasePressures];
       
end


end

